El concepto de entropía
La entropía se puede interpretar como un indicativo de la complejidad en una serie de datos aleatorios \(\{X_1, X_2, \ldots, X_N\}\). La entropía usualmente se calcula utilizando una función de masa de probabilidad \(p_j\) que cumple con las siguientes propiedades:
- \(p_j \ge 0, \, j \in 0,1, \ldots, N\)
- \(\sum_{j=1}^N{p_j} = 1\)
R tiene múltiples paquetes y funcionalidades que permiten estimar la pmf de un conjunto de datos como el descrito anteriormente. El histograma es una herramienta que permite estimar la pmf de un conjunto de datos. En el ejemplo que sigue se muestra la forma de estimar la pmf en R:
set.seed(1234) # Para hacer el análisis reproducible
datos <- rnorm(512,0,1) # Se generan 512 valores normales
histogram <- hist(datos, plot=FALSE) # Se calcula el histograma
pmf <- histogram$counts/sum(histogram$counts) # Se calcula la pmf
sum(pmf) # Se verifica que cumpla con las propiedades
[1] 1
Ejercicios
- Estimar la
pmfutilizando el utilizando los paquetesASHyKernSmoooth.
library(ash)
a <- ash1(bin1(datos, ab = c(min(datos), max(datos))))
[1] "ash estimate nonzero outside interval ab"
plot(a, type = "l")
sum(diff(a$x) * (head(a$y, -1) + tail(a$y, -1)))/2
[1] 0.9973948
library(KernSmooth)
k <- bkde(x = datos, range.x = c(min(datos), max(datos)))
plot(k, type = 'l')
- ¿Cuál es la ventaja de utilizar los métodos anteriores sobre el histograma?
El histograma es el método más viejo y menos sofisticado para estimar la densidad. ASH utiliza el método de average shifted histogram, el cual es más suave que el histograma y evita la sensibilidad a la elección del origen, pero sigue siendo computacionalmente eficiente.
KernSmooth utiliza el método de kernel density estimation. El kernel es una función simétrica, generalmente positiva, que se integra en 1. El enfoque de kernel density estimation supera la discreción de los enfoques del histograma al centrar una función de kernel suave en cada punto de datos y luego sumar para obtener una estimación de densidad.
- Utilizando el comando
histy los paquetesASHyKernSmoothverifique el tiempo requerido para estimar la densidad de una serie de datos Gaussianos con \(\mu=1\) y varianza \(\sigma^2=1\) y longitudes \(N=2^i, \, i=8,9,10, 11, \ldots 16.\) (es necesario incluir un gráfico enhighcharter)
dataSeries <- list()
for (i in 8:16){
dist <- rnorm(2^i, mean = 1, sd = 1)
dataSeries[[(i-7)]] <- dist
}
dataAll = do.call(rbind, dataSeries)
tiempoHist <- list(NA, 9)
for (i in 1:9){
start <- Sys.time()
hist(dataAll[i,], plot=FALSE)
end <- Sys.time()
runtime <- end - start
tiempoHist[[i]] <- runtime
}
th <- do.call(rbind, tiempoHist)
nombre <- rep("Hist", 9)
n_i <- (8:16)
th <- cbind(th, nombre, n_i)
tiempoASH <- list(NA, 9)
for (i in 1:9){
start <- Sys.time()
ash1(bin1(dataAll[i,], ab = c(min(dataAll[i,]), max(dataAll[i,]))))
end <- Sys.time()
runtime <- end - start
tiempoASH[[i]] <- runtime
}
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
tash <- do.call(rbind, tiempoASH)
nombre <- rep("ASH", 9)
n_i <- (8:16)
tash <- cbind(tash, nombre, n_i)
tiempoKS <- list(NA, 9)
for (i in 1:9){
start <- Sys.time()
bkde(x = dataAll[i,], range.x = c(min(dataAll[i,]), max(dataAll[i,])))
end <- Sys.time()
runtime <- end - start
tiempoKS[[i]] <- runtime
}
tks <- do.call(rbind, tiempoKS)
nombre <- rep("KS", 9)
n_i <- (8:16)
tks <- cbind(tks, nombre, n_i)
tiempos <- rbind(th, tash, tks)
tiempos.df <- as.data.frame(tiempos)
names(tiempos.df)[1] <- "time"
names(tiempos.df)[2] <- "method"
names(tiempos.df)[3] <- "i"
tiempos.df$time <- as.numeric(tiempos.df$time)
head(tiempos.df)
time method i
1 0.01270819 Hist 8
2 0.01890492 Hist 9
3 0.02113605 Hist 10
4 0.02066207 Hist 11
5 0.01725292 Hist 12
6 0.01745105 Hist 13
library(highcharter)
Gráfica de tiempos
hchart(tiempos.df, "column", hcaes(x = i, y = time, group = method))
La entropía utilizando el histograma
Volviendo de nuevo al ejemplo anterior, podemos estimar la entropía de Shannon, utilizando la pmf obtenida mediante el histograma y así obtener un estimador empírico de la entropía de Shannon. A continuación mostramos la forma de obtener la entropía de una serie de datos obtenida en ventanas independientes o contiguas de longitud \(512\):
set.seed(1234)
datos <- rnorm(32768)
wLength <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")
noVentanas <- length(datos)/wLength
entropies <- numeric(noVentanas)
index <- numeric(noVentanas)
for(i in 1:noVentanas)
{
dataW <- datos[wLength*(i-1)+1:wLength*i]
histo <- hist(dataW, breaks=8,plot=FALSE)
pmf <- histo$counts/sum(histo$counts)
entropies[i] <- -1*sum(pmf*log(pmf))
index[i] <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")
Preguntas
- ¿Porqué existen valores discontinuos en la entropía?
Porque en algunos de los resultados de la pmf el resultado es 0 y al calcular la entropía obtenemos valores NaN.
- ¿Con qué código soluciona el problema de las discontinuidades?
Podemos excluir los valores NaN con el parámetro na.rm
entropia <- function(datos, wLength)
{
noVentanas <- length(datos)/wLength
entropies <- numeric(noVentanas)
index <- numeric(noVentanas)
for(i in 1:noVentanas)
{
dataW <- datos[wLength*(i-1)+1:wLength*i]
histo <- hist(dataW, breaks=8,plot=FALSE)
pmf <- histo$counts/sum(histo$counts)
entropies[i] <- -1*sum(pmf*log(pmf), na.rm = TRUE)
index[i] <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")
}
entropia(datos, wLength)
- Ahora calcule la entropía de una serie de datos normales, denotados por \(X_t\) (con \(\mu=0\) y \(\sigma=2\)) pero ahora añadanle (súmenle) una segunda función \(r_t\), es decir, hallen la entropía de la serie \(Y_t = X_t+r_t\) con \(r_t\) definida por: \[ r_t = \begin{cases} \sigma/4 & t\ge 16384\\ 0 & \mbox{otro caso} \end{cases} \]
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/4, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
entropia(datos, wLength)
- Repitan el paso \(3\) para \(r_t\) dada por: \[ r_t = \begin{cases} \sigma/2 & t\ge 16384\\ 0 & \mbox{otro caso} \end{cases} \]
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/2, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
entropia(datos, wLength)
- Repitan el paso \(3\) para \(r_t\) dada por: \[ r_t = \begin{cases} \sigma & t\ge 16384\\ 0 & \mbox{otro caso} \end{cases} \]
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
entropia(datos, wLength)
- ¿Tiene algún efecto la longitud del salto en \(r_t\) en la forma de la entropía? Explique.
Entre mayor sea el valor de sigma se presentan menos picos en la entropía.
- ¿Qué sucede ahora si \(r_t\) es de la forma: \[ r_t = \begin{cases} \sigma & 15872 \le t\le 16896\\ 0 & \mbox{otro caso} \end{cases} \]?
r_t <- function(t, sigma) ifelse(15872 <= t & t<= 16896, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength <- 512
entropia(datos, wLength)
Los picos de los valores intermedios disminuyen.
- Repita los pasos 3-7 pero ahora usando la entropía de Harvda con parámetro \(\alpha=3\) y \(\alpha=9\). ¿Qué efecto tiene \(\alpha\)?
havrda <- function(pmf, alpha) sum(pmf^alpha - 1) / (2^(1 - alpha) - 1)
entropiaHavrda <- function(datos, wLength, alphaVal)
{
noVentanas <- length(datos)/wLength
entropies <- numeric(noVentanas)
index <- numeric(noVentanas)
for(i in 1:noVentanas)
{
dataW <- datos[wLength*(i-1)+1:wLength*i]
histo <- hist(dataW, breaks=8,plot=FALSE)
pmf <- histo$counts/sum(histo$counts)
entropies[i] <- havrda(pmf, alpha = alphaVal)
index[i] <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", xlab="Tiempo, t", ylab="Valores de entropía")
}
- Paso 3
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/4, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
\(\alpha = 3\)
entropiaHavrda(datos, wLength, 3)
\(\alpha = 9\)
entropiaHavrda(datos, wLength, 9)
- Paso 4
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/2, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
\(\alpha = 3\)
entropiaHavrda(datos, wLength, 3)
\(\alpha = 9\)
entropiaHavrda(datos, wLength, 9)
- Paso 5
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
\(\alpha = 3\)
entropiaHavrda(datos, wLength, 3)
\(\alpha = 9\)
entropiaHavrda(datos, wLength, 9)
- Paso 7
r_t <- function(t, sigma) ifelse(15872 <= t & t<= 16896, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
\(\alpha = 3\)
entropiaHavrda(datos, wLength, 3)
\(\alpha = 9\)
entropiaHavrda(datos, wLength, 9)
Podemos observar que entre menor sea el valor de \(\alpha\) el valor de la entropía disminuirá.
Entropía utilizando ventanas deslizantes
El cálculo de la entropía por ventanas independientes dada arriba resulta útil en casos en dónde la función no tiene dependencia en los valores futuros (descorrelacionadas). Para el caso de funciones correlacionadas, el cálculo de la entropía por ventanas deslizantes traslapadas resulta útil para descubrir ciertas fenomenologías en los datos. El cálculo por ventanas deslizantes de tamaño \(W\) se realiza sobre una secuencia de datos \(X_1, X_2, \ldots, X_N\). La ventana (\(W\le N\)) se va deslizando sobre los datos con factor \(\Delta\) y de esta forma subconjuntos de los datos \(X_i\) toman la siguiente forma: \[ X(m; W, \Delta) = x_j \times \Pi(\frac{t-m\Delta}{W}-\frac{1}{2}), \] donde \(m\Delta \le j \le m\Delta + W\) y \(m=0,1,2, \ldots\). Finalmente se puede graficar \(nW + \Delta, n=1,2,3, \ldots\) contra las entropías y verificar algún patrón en los datos.
Ejercicios
Implementar en
Rla metodología del cálculo de la entropía de Harvda normalizada por ventanas deslizantes. La función debe tener la formaharvda_deslizante(datos, w.length=512, s.factor=10, a.parameter=0.8, ent.type=c("hist", "ash", "kern")), dondew.lengthes la longitud de la ventana,s.factores el factor de deslizamiento ya.parameteres el parámetro \(\alpha\) de la entropía de Harvda. Además la función puede calcular la entropia usando el histograma, por el método ash o por alguna metodología kernel (con el parámetroent.type).Aplicar la entropía calculada con los datos generados anteriormente, es decir:
Los pasos \(3,4,5\) y \(7\) en donde los datos se dan como \(X_t+r_t\).
havrdaNormalizada <- function(pmf, N, a.parameter)
{
sum(pmf^a.parameter - 1) / (N^(1 - a.parameter) - 1)
}
havrdaDeslizante <- function(datos, w.length = 512, s.factor = 10, a.parameter = 0.8, ent.type=c("hist", "ash", "kern"))
{
N <- length(datos)
m <- 0
n <- 1
inf <- 0
sup <- 0
entropies <- c()
index <- c()
inf <- m * s.factor + 1
sup <- m * s.factor + w.length
while (sup <= N)
{
dataW <- datos[inf:sup]
pmf <- switch (ent.type,
"hist" = {hist <- hist(dataW, breaks = 1000, plot = FALSE)
hist$counts / sum(hist$counts)},
"ash" = {invisible(capture.output(hist <- ash1(bin1(dataW))))
hist$y / sum(hist$y)},
"kern" = {hist <- bkde(dataW)
hist$y / sum(hist$y)})
entropies <- append(entropies, havrdaNormalizada(pmf, N, a.parameter))
index <- append(index, n*w.length + s.factor)
m <- m + 1
n <- n + 1
inf <- m * s.factor + 1
sup <- m * s.factor + w.length
}
return(list(index, entropies))
}
- Paso 3
library(zeallot)
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/4, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "hist")
plot(index, entropies, type = "l", main="Usando hist", xlab="Tiempo, t", ylab="entropy")
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "ash")
plot(index, entropies, type = "l", main="Usando ash", xlab="Tiempo, t", ylab="entropy")
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "kern")
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")
- Paso 4
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/2, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "hist")
plot(index, entropies, type = "l", main="Usando hist", xlab="Tiempo, t", ylab="entropy")
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "ash")
plot(index, entropies, type = "l", main="Usando ash", xlab="Tiempo, t", ylab="entropy")
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "kern")
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")
- Paso 5
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "hist")
plot(index, entropies, type = "l", main="Usando hist", xlab="Tiempo, t", ylab="entropy")
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "ash")
plot(index, entropies, type = "l", main="Usando ash", xlab="Tiempo, t", ylab="entropy")
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "kern")
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")
- Paso 7
r_t <- function(t, sigma) ifelse(15872 <= t & t<= 16896, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "hist")
plot(index, entropies, type = "l", main="Usando hist", xlab="Tiempo, t", ylab="entropy")
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "ash")
plot(index, entropies, type = "l", main="Usando ash", xlab="Tiempo, t", ylab="entropy")
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "kern")
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")
- ¿Tiene algún efecto el parámetro \(a\) en la forma de la entropía? ¿Tiene alguna ventaja calcular la entropía de Harvda por otro método diferente al histograma?
Comparamos utilizando kern y modificando el valor de \(\alpha\)
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
\(\alpha\) = 1.5
c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 1.5, ent.type = "kern")
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")
\(\alpha\) = 3
c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 3, ent.type = "kern")
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")
\(\alpha\) = 5
c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 5, ent.type = "kern")
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")
\(\alpha\) = 7
c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 7, ent.type = "kern")
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")
\(\alpha\) = 9
c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 9, ent.type = "kern")
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")
\(\alpha\) traslada la función sobre el eje y también disminuye el tamaño de la función.
- Además del histograma y los estimadores tipo kernel, existen otros métodos para estimar la distribución de una serie de datos. Investigue: ¿en qué consiste la entropía de permutación?
Según Judge y Henry (2019) la entropía de permutación es una herramienta robusta de series de tiempo que proporciona una medida de cuantificación de la complejidad de un sistema dinámico al capturar las relaciones de orden entre los valores de una serie de tiempo y extraer una distribución de probabilidad de los patrones ordinales.
Referencias
Deng, H. & Wickham, H.. (2011). Density estimation in R.
Henry, M. (2020). Permutation Entropy. Aptech
Barnier, J. (2020). Rmdformats: HTML Output Formats and Templates for Rmarkdown Documents. https://github.com/juba/rmdformats.